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ABSTRACT 


An explicit numerical solution of the compressible Navier-Stokes 
equations is applied to the thermodynamic analysis of supercritical 
oxygen in the Apollo cryogenic storage system. The wave character is 
retained in the conservation equations which are written in the basic 
fluid variables for a two-dimensional Cartesian coordinate system. 
Control-volume cells are employed to simplify imposition of- boundary 
conditions and to ensure strict observance of local and global conser- 
vation principles. Non-linear real-gas thermodynamic properties respon 
sible for the pressure collapse phenomonon in supercritical fluids are 
represented by tabular and empirical functions relating pressure and 
temperature to density and internal energy. Wall boundary conditions 
are adjusted at one cell face to emit a prescribed mass flowrate. Elec 
trical heater input is treated as localized internal heat generation, 
a fraction of which may be radiated to the walls where it is added to 
the prescribed boundary heat flux. The effect of "tank stretch" on 
dP/dt is included as out-of-plane fluid expansion. Scaling principles 
are invoked to achieve acceptable computer execution times for very low 
Mach number convection problems. Detailed simulations of thermal 
stratification and fluid mixing occurring under low acceleration in the 
Apollo 12 supercritical oxygen tank are presented which model the 
pressure decay associated with de-stratification induced by an ordinary 
vehicle maneuver and heater cycle operation. 
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INTRODUCTION 


Origin of problem 

Gaseous oxygen required for life-support and electrical power 
generation aboard the Apollo Command and Service Module is stored in 
super-insulated double-wall spherical tanks which are a portion of the 
storage and supply system. In these tanks, oxygen is maintained at 
cryogenic temperatures and at pressures somewhat above the critical 
pressure. (See Table 1 for system operating parameters.) Under these 
conditions, the oxygen is in the super-critical state which results in 
a high-density single-phase compressible fluid suitable for expulsion 
under zero-gravity conditions. 


Storage tank pressure is maintained relatively constant during 
fluid withdrawal by increasing the fluid temperature using an internal 
filament-type electrical heater. Prior to and including Apollo 13, the 
energy supplied by this heater was dispersed throughout the fluid by 
means of internal mixing fans. However, following the failure of the 
oxygen system during the Apollo 13 mission, the mixing fans were re- 
moved from subsequent oxygen tanks to minimize the possible combustion 
hazard. 

[21 

Numerical studies by Kamat and Abraham have shown that a 
substantial decrease in pressure can result from fluid mixing if the 
heater input is not well distributed beforehand. The decrease in pres- 
sure can be of sufficient magnitude to cause the fluid to return to a 
two-phase condition. Since natural convection would have to be relied 
upon to limit the temperature concentrations, a considerable effort was 
initiated to refine the understanding of the thermodynamic and fluid 
dynamic behavior of locally-heated supercritical oxygen stored in a 
low acceleration environment. 

Stratification and Pressure Collapse 

Due to the low thermal conductivity of supercritical oxygen, 
electricial heaters generate local heat concentrations. Under low 
accelerations, these concentrations do not disperse rapidly. Without 
a mechanical means of mixing the heated fluid with the surrounding 
cold fluid, subsequent heater cycles result in increasingly severe 
heat concentrations, or "thermal stratification." 

Thermophysical properties of oxygen in the vicinity of the criti- 
cal point are strongly nonlinear and therefore deviate from the pro- 
perties of an ideal gas. One of the ramifications of the specific heat 
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nonlinearity is that the bulk fluid temperature is not equal to the 
equilibrium temperature which would exist if the fluid were mixed or 
otherwise brought to thermal equilibrium. 

The equation of state is also nonlinear in this region and is 
sensitive to small changes in temperature. The result is that a 
modest drop in the bulk temperature brought about by fluid mixing from 
a stratified condition can cause a substantial decrease in pressure or 
"pressure collapse". Kamat and Abraham have shown that the pres- 
sure decreases monotonically with mixing so that the final pressure is 
always below any intermediate pressure. It is possible for the 
collapse pressure to drop below the critical pressure in which case 
the fluid becomes sub-critical, a condition in which liquid and gaseous 
phases may co-exist. This condition must be avoided if a uniform single 
phase fluid expulsion is to be maintained. 


This paper describes an explicit numerical solution of the com- 
pressible conservation equations that govern natural convection in 
compressible fluids. Real fluid thermodynamic relations are employed 
so that the pressure collapse phenomenon can be observed in super-criti- 
cal fluids. A number of additional effects have been incorporated in 
the solution for engineering application including "tank stretch," 
heater thermal mass, and heater radiation. 


Comparisons are made between the results of the present solution 
and Apollo 12 flight data occurring at 8 hours GET. Several other 
verification test cases also are presented to demonstrate the program 
capability. 

Table 1 Cryogenic System Operational Parameters - Oxygen 
(Apollo 14 And Subsequent) 


Stored Fluid Weight (100% indication) 
Usable Fluid Weight 
Tank Volume 

Normal Operating Pressure 

Pressure Switch Deadband (Min) 

Total Heater Power 

(prior to Apollo 14) 

Bulk Fluid Temperature 

4r 

Data taken from Reference 1 


330.1 lb 

323.5 lb 

4.75 ft 3 

900 + 35 psi 

30 psi 

434.8 B/hr 
(528.6) 

160 to 530 °R 
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APPROACH 


Under low acceleration conditions, the potential pressure collapse 
depends primarily upon the amount of heat concentrated in the portion 
of fluid surrounding the heater. The heater cycle time, however, de- 
pends upon the actual pressure rise and decay rates which, in turn, 
are related to the processes taking place in the heater thermal boun- 
dary layer and to other effects such as tank stretch. At 95% of tank 
quantity, tank stretch decreases the magnitude of the pressure rate and 
therefore the cycle time by almost a factor of two^K 

A detailed knowledge of the temperature, density, and fluid 
velocity distributions is necessary in order to describe the heating 
and mixing processes and the resulting thermodynamic behavior of super- 
critical oxygen in the Apollo cryogenic storage tanks. This informa- 
tion can be obtained only by simultaneously solving the conservation 
equations of mass, momentum, and energy. These equations are coupled 
nonlinear partial differential equations which are not amenable to 
analytical solutions for any but the most restrictive of problems. 
However, finite-difference numerical techniques have been developed 
which allow these equations to be solved in their entirety. 

The general solution procedure is to integrate numerically the 
conservation equations from a given set of initial conditions subject 
to the various boundary conditions imposed by the physical system. 
Although the fluid motion in the Apollo storage tanks is three-dimen- 
sional, a two-dimensional solution of the conservation equations is 
adequate to resolve the basic mechanisms which produce stratification 
and mixing. 

The basic model consists of a two-dimensional fluid slab of unit 
depth to which the conservation equations are applied. The momentum 
equations include the gravitational body force terms so that natural 
convection may develop. The basic provisions necessary to model the 
Apollo storage tank are: two acceleration components, a localized 

internal heat source, external heat leak, and a fluid outlet port. 
Non-linear thermodynamics require that real-gas properties be used. 

To adequately describe the actual pressure rise rate for engineering 
purposes, a number of refinements to the idealized model were required 
which include accounting for the thermal mass of the heater, heater 
radiation to the tank wall, and tank stretch. 

Acceleration Components 

In general, two acceleration components are required to describe 
the body forces acting on the fluid. Stratification develops under a 
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uniform low-level acceleration in a constant relative direction. Such 

an acceleration arises from the centripetal acceleration (w x oj x r) 
associated with the three revolutions per hour roll rate used for pas- 
sive thermal control (PTC) . 

Maximum mixing resulting from a G-spike, (e.g., RCS thruster fir- 
ing) occurs when the direction of the G-spike is normal to the density 

-+■ -y 

gradient. A tangential acceleration (w x r) results from a change in 
the roll rate during spin-up to or de-spin from the PTC mode. The 
effects of coriolis acceleration, non-constant radius vector, and 
planetary gravity gradient are insignificant compared to the accelera- 
tion forces discussed above and are not considered in this analysis. 

Internal Heater 


The heater element is simulated as an internal energy source at 
an appropriate location in the fluid volume. No velocity boundary 
conditions are attempted at these "heater cells" due to the computa- 
tional difficulty in resolving the flow boundary layer. The heater 
does not represent the heater cylinder itself but rather the heated 
fluid sheath surrounding the heater cylinder. The number of heater 
cells is selected such that the combined volume is equal to the 
"effective boundary layer volume". This effective volume has been 
postulated in order to explain observed pressure rise rates, since the 
grid spacing of the numerical solution is too coarse to resolve the 
actual thermal boundary layer. The effective boundary layer volume 
was derived on the basis of empirical flight data. 

Heater Thermal Mass Effect 


Heater-on operation is modeled as a prescribed heat generation 
rate within the heater cells . When steady-state heating conditions are 
approached, a power balance at the heater implies that the impressed 
electrical power is converted directly to a heating rate of the fluid 
sheath surrounding the heater. The heater thermal mass absorbs heater 
power until heat conduction into the fluid reaches steady-state. For 
the present purposes, this build-up rate is approximated as an equi- 
valent linear ramp up to the steady-state rate which occurs during 


the interval T, 


(See Nomenclature.) Defining C, 


L • V JCC I'lUlUClit.iaLUtCi / 1/U1.XIIXU5 

lag lag 

as a ramp function, the internal heat generation rate then is given 

by: . . 

Q = C, Q, 

lag heater 


( 1 ) 
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After the heater is turned off, the heater cylinder rejects the 
heat stored in its thermal mass. Assuming that the heater temperature 
returns to the previous low temperature, the inverse linear ramp was 
used for flux decay. The internal heat generation rate during the 
heater-off ramp then becomes 


Q = (1 - C, ) (1 

lag heater 


( 2 ) 


Heater Radiation 


During the latter portions of the Apollo mission, the heater 
temperature rises high enough that thermal radiation plays an import- 
ant heat transfer role. At a specified time in the mission, the energy 
radiated from the heater may be represented as a constant fraction of 
the heater power not absorbed by the thermal mass: 


Q , = C . (C, (1 ) 

rad rad lag heater 

and the rate of heat entering the fluid becomes: 

Q = (1 - C ) C Q 

rad lag heater 


Outlet Port 


(3) 

(4) 


A fluid outlet port is modeled at the periphery of the fluid 
volume by specifying a fluid velocity at a location on the boundary 
such that the prescribed mass withdrawal rate occurs. This same 
velocity is used to convect momentum and energy from the system. 

It was necessary to relocate the outlet port from its actual 
position to a position lying in the plane of the model. The flow 
distortion introduced by this relocation appears to be negligible for 
the flowrates presently being considered. 


Heat Leak 


Heat entering the fluid after passing through the super-insulation 
surrounding the storage tank is called heat leak and for the Apollo oxy- 
gen tank has a nominal value of about 25 BTU/hr. The quantity enter- 
ing the fluid model is proportioned according to the volume ratio C . 

Heat leak is imposed as a uniformly distributed heat flux over the 
exposed fluid boundary. The decision to specify the heat flux at the 
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boundary as opposed to specifying the boundary temperature was based 
upon the fact that the rate of heat leak into the fluid is more 
accurately known than the wall temperature of the inner tank. Global 
conservation principles also are more easily satisfied. Since oxygen 
is primarily transparent in the region, the energy radiated from the 
heater cylinder is not absorbed until it reaches the tank wall from 
which it enters the fluid by conduction. 

Tank Stretch/Line Compression 

The elasticity of the thin-wall tank permits a volume expansion 

-5 3 

with pressure (dV/dP) of about 3.6 x 10 ft /psi at the normal 
operating pressure. At high tank quantities (80%-100%) fluid pressure 
is very sensitive to density, and variations in tank volume have a 
significant effect on the pressure rise rate p3 - 16] decreasing 
the magnitude of dP/dt by over 40% at 95% quantity. Tank stretch is 
modeled as an out-of-plane fluid expansion by permitting the unit z 
dimension to vary with pressure according to dV/dP. 

[3] 

Seto derived an expression for the effect on dP/dt of fluid 
compression in the external plumbing. Tests with the equilibrium 
tank model described in this reference show that this effect is of the 
same order as the effect of tank stretch. For the present purpose, 
line compression is represented as an increase in the tank stretch 
factor dV/dP. 
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GOVERNING EQUATIONS 


General Discussion 


The equations governing the conservation of mass, momentum, and 
energy in a fluid system may be formulated in either Eulerian or 
Lagrangian coordinates. The Lagrangian formulation fixes the coordinate 
system in the fluid mass and is convenient for problems involving a 
free surface. However, accuracy is seriously impared with time as the 

coordinates deform and move with fluid convection . In the Eulerian 
formulation, the coordinate system is fixed in the fluid volume so that 
the fluid moves through the coordinate system. For this reason, the 
Eulerian formulation is generally superior for complex convection 
problems and is used for the present analysis. 


[A] 

Richtmyer and Morton and others recommend the use of conser- 
vation equations in the "conservative" or "divergence" form which 
preserves the conservation principles when solved at a finite number of 
discrete points. If the conservative form is not used, computational 
sources and sinks can appear as fluid is convected from one cell to 
the next, their origin being the non-constant coefficients appearing in 
front of the derivatives. For example, the mass leaving a cell through 
one of its faces does not necessarily appear in total in the adjacent 
cell. Although these errors are small, they occur at each cell inter- 
face at each time step and can accumulate with time. The propagation 
of these errors eventually can lead to computational instability. Com- 
putationally non-conservative equations are derived by performing flux 


balances on an infinitesimal control volume . These equations can be 
converted to conservative form by adding the continuity equation to 
each of the other three conservation equations. The divergence form is 
developed directly from surface and volume integrals of flux vectors 

employing the Gauss Divergence Theorem^. 


A very clean form of the general conservation equations was ob- 
tained with minor modification from Goodrich These equations are 

in conservative form and are vrritten in terms of the basic fluid vari- 
ables for a two-dimensional Cartesian coordinate system. The non- 
dimensional form was not used because of the difficulty in assigning 
characteristic reference values for non-linear real fluid thermodynamic 
properties . 
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The general equations governing a compressible viscous fluid with body 
forces and internal heat source are: 


Mass Continuity 


dp , d (pu) , 3 (pv) 

3t 3x 3y 


X-Component Momentum 


3 (pu) 
3t 


— (puu + P + T ) 
3x xx 


3y 


(puv 


+ T ) 
xy 


Pg, 


Y-Component Momentum 


3 (pv) 

3t 


+ T — (pvu + T ) 
3x K xy 


+ — (pvv + P + T ) 

3y yy 


pg. 


Energy 


3 (pE) 3 r 

3 1 + 3^ 

+ hfi 


where: 
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E = e + ^(u^ + v^) 

P = P(p,e) 

T = T(p,e) 

, 3T 

q x = " k 3^ 

. 3T 

q = -k -r— 

y 3y 

2 ,„3u 3v, 

t = - -ry (2r — ) 

xx 3 3x 3y 

2 ,_3v 3u, 

T = - -rp (2- — ) 

yy 3 3y 3x 


,3u 3v, 

T xy - - V( 87 + 


(5) 


( 6 ) 


(7) 


( 8 ) 

(9) 

( 10 ) 
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( 12 ) 
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The form of the energy equation in which temperature is the inte- 
grated variable ’ ^ is not appropriate for supercritical fluid 
analysis since this form assumes a constant specific heat (C ). 


Program Formulation 

The above equations are more general than necessary for the present 
analysis so that a number of non-critical terms were eliminated to 
improve computational speed. 

An advantage in using the internal energy form of the energy equation 
is that the fluid temperature is almost entirely dependent upon the 
specific internal energy and only slightly dependent upon density. 
Therefore, over a relatively wide range of pressures, the relation T= 

T(p ,e) was replaced by the computationally simpler function involving 
only a one-dimensional table interpolation: 

T = T (e) (14) 

roi 

Figure 1 shows this function (obtained from Weber's 1 J data) 
platted at 55, 60 and 65 atmospheres to indicate the error incurred by 
neglecting pressure variation in the vicinity of the critical point. 

The form of the equation of state shown by eqn. (10) was not 
available, so that it was necessary to use the equivalent (though for 
this application computationally less desirable) relation 

P = P (P, T) (15) 

T91 

which was available in Stewart's equation of state . The relations 
selected for pressure and temperature are acceptable in that an itera- 
tive procedure is not required. However, if pressure is to be computed 
from temperature, the above approximation to the temperature-energy 
function may be more critical than if temperature is used only for heat 
conduction. 

The general governing equations assume only that Stokes hypothesis 
regarding the viscosity coefficients holds, and that thermal radiation/ 
absorbtion effects are insignificant. For the present application, a 
number of simplifications to these equations are made. First, the 
velocities developed in the low acceleration environment are so small 
that the kinetic energy terms in the energy equation may be neglected. 
Second, the low velocities coupled with the low viscosity of oxygen 
make the viscous terms in the energy equation negligible. For compu- 
tational simplicity thermal conductivity and viscosity are assumed 
constant. The remaining viscous terms are further simplified by 
assuming that the fluid is incompressible as far as viscous dissipation 
is concerned. 
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Incorporating these simplifications in the general equations and 
rearranging, the governing equations used in this analysis are written. 


Continuity 


3p _ 3 (pu) 3 (pv) 

3t 3x 3y 


X-Component Momentum 

3 (pu) _ 3 (puu) _ 3 (puv) 

3t 3x 3y 

Y-Component Momentum 

3 (pv) _ 3 (pvu) _ 3 (pvv) 

3t 3x 3y 

Energy 

3 (pe) _ _ 3 (pe-fP ) 3 (pe+P) 
3t 3x 3y 

Thermodynamic Relation 

T = T(e) 

Equation of State 

P = P(P,T) 


ajp 

3x 


+ Pi 


3P 

3y 


+ pg y 



+ y (■ 


3^ 

3x" 


3y 


.3^v 3^ 

+ ' > ^ + -y 



Typical initial conditions are: 


u(x,y,o) = 0 
v(x,y,o) = 0 
P(x,y,o) = P o (x,y) 
T(x,y,o) = T q 


(specified to balance body forces) 


(16) 

£) ( 17 ) 

£> < 18 > 

(19) 

( 20 ) 
( 21 ) 


( 22 ) 



Boundary conditions for a closed tank take the form: 
u(w,y,t) = u(x,w,t) = 0 

v(w,y , t) = v(x,w,t) = 0 (23 

|p (w,y , t) = |p- (x,w,t) = - (heat leak) 

where w indicates a value of x or y at any wall. 

The addition of boundary conditions defining a heater and an out- 
let port is discussed more thoroughly in a following section. Basic- 
ally, however, the heater is a region of fluid in which internal heat 
generation is specified, and the outlet port is defined by a normal 
velocity at a section of the tank wall such that a prescribed mass 
withdrawal rate occurs. This exit velocity convects mass, momentum, 
and energy from the system. 
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FINITE DIFFERENCE FORMULATION 


Discussion of Numerical Techniques 

A wide variety of finite-difference schemes are available for 
evaluating the governing equations, and Richtmyer and Morton ^ pro- 
vides a good source reference. The AIAA reprint series contains 

an extensive bibliography and a collection of the more interesting 
recent papers relating primarily to high-speed flow. Numerical methods 

discussed in Ralston and Wilfe^^, and Cheng^^] summarizes the funda- 
mental principles relating to the numerical solution of the Navier-Stokes 
equations . 

The numerical integration of the governing equations is performed 
at a finite number of discrete points located throughout the fluid 
volume. The difference equations are obtained by replacing the partial 
derivatives with suitable finite-difference approximations typically 
derived by Taylor-series expansions in space and time. Alternately, 
the difference equations can be derived directly from fundamental con- 
servation principles applied to a fluid control volume ^ 12 \ This 
method avoids taking the limit as AV-K) to form the differential equa- 
tions followed by the reverse process of discretization. The two 
methods are basically equivalent; however, the latter is quite useful 
for visualizing and formulating conservative- form difference equations 
particularly in curvilinear coordinate systems and parameter spaces. 

The time variable is discretized and is given by t n = n At. The 
fluid state is advanced from the (n) time plane to the (n+1) time plane 
by integration. When the time derivative is approximated by a forward 
difference, all information used to advance the state to (n+1) At is 
available at the n^-* 1 time plane. This is the explicit formulation and 
is equivalent to a step integration procedure. The two-step Lax-Wen- 
droff sc heme generates provisional values at (n + h) At which are 
used to advance to the (n+1) time plane. This procedure centers the 
time difference which gives second order accuracy in time. 

Various other numerical schemes are available which use weighted 
averages of data obtained from several time planes (n-1) , (n) , and 
( n +l) to obtain the updated value at (n+1) . Schemes which require data 
at the (n+1) time plane to advance the time to (n+1) At are implicit 
and require the simultaneous solution of high-order systems of equations. 
Efficient relaxation techniques make the implicit formulation quite 
powerful for certain types of problems. 
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Most implicit schemes are unconditionally stable for any time step, 
whereas, explicit schemes are stable only when rather stringent stability 
conditions are observed. However, explicit schemes are straight forward 
and require the least computation per time step. The advantage of one 
scheme over the other depends upon the rate at which the fluid properties 
are changing. Explicit schemes are preferred for time-dependent prob- 
lems in high-speed flow and for problems in which wave propagation is 
important. Implicit schemes are effective for problems in low speed 
flow and certain steady-state problems. The use of scaling principles 
to transform the low speed flow problem into one in which the transport 
mechanisms occur very much faster than in real time appears to be effec- 
tive in broadening the class of problems which may be solved efficiently 
using the explicit formulation. 

Problems can arise in achieving stability in the explicity-formulated 

[■41 

conservation equations. Lax and Wendroff 1 used second-order accurate 
centered time differences. Rusanov ^ ^ ] added numerical damping terms to 

the space differences. Goodrich has shown that the numerical damp- 
ing terms of Rusanov can be represented by a weighted biasing of the 
convective difference terms; a method which is similar to the upstream 

differencing technique used below. Richtmyer^^ has suggested adding 

pseudo-viscosity terms of the form a | | to stabilize the momentum 

equations . 


Technique Employed 

The technique employed in this analysis is patterned after the method 

[4 p 292] 

of Courant, Fredricks, and Lewy (1929) for the wave equation ’ p 
In this method, the momentum equations are advanced to the (n+1) time 

plane first, and the updated velocities u n+ \ v n+ ^ are then used to 
advance the continuity and energy equations. By performing the integra- 
tion in two steps, the continuity and energy equations appear to be semi- 
implicit, although they are effectively explicit since the necessary 
information is available from previous calculation. However, when this 
method is used with centered space differences, physically unrealistic 
temperature profiles result as fluid convection takes place. Richtmyer 

[4, p 292] cre( j^ ts the cure f or this problem to Lelevier. This solution 
replaces the centered space differences used for the convective terms 
by forward or backward space differences as the sign of the convecting 
velocity is negative or positive. This procedure is quite common and goes 
by a number of names including upwind or upstream differencing and 
donor-cell differencing. These differences, however, are only of first 
order accuracy. 
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As applied to the grid system used in this analysis, upstream 
differencing has the effect of defining the values of the convected 
properties at the control- volume interface as being the values of the 
upstream cell. Central differencing implies averages of the adjacent 
cell values at the interface. 

Grid System 

Grid points are uniformly spaced in the x and y directions at 
constant intervals Ax = Ay = L. It is convenient to visualize an 
elemental control volume AV = Ax Ay Az = L , referred to as a cell, 
surrounding each grid point. According to Cheng , "The fluxes must 

be evaluated on the cell boundary while the conserved quantities are 
determined only as averages over the cell." 

The two- dimensional fluid volume is characterized by several 
hundred of these cubic cells arranged in a plane. The circular con- 
figuration of the oxygen tank crosssection is approximated in a step- 
wise fashion by removing cells from comers of the rectangular cell 
grid. For example see Figure 9. 

Individual cells are identified by the indicies (i,j) in the x- 
and y- direction as shown in Figure 2 so that the physical position 
of the center of cell (i,j) is located at 

(x^y^) = (i-^Ax, (j-Js)Ay (24) 

The integrated fluid properties p, pu, pv, and pe are identified 

with each cell center and represent the average properties over the cell 
volume. The value of a property at the center of cell (i,j) is design- 
ated for example as 





(25) 


For mathematical consistency, the velocities and the thermodynamic 

properties of the fluid in the cell also are defined at the cell center 
so that the following relationships can be used: 


u 


ij 




ll»7 



(26) 
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Pv 



ij 


T. . = T 
ij 


(pe ij /p ij ) 



(p u,V 


It is convenient to define properties, certain gradients, and 
convection velocities at the cell walls which lie midway between grid 
points, and which are indicated by half-subscripts. Convection velo- 
cities at these faces are evaluated by linear interpolation: 


u. k a = ^(u- , . + u ) 
i“ 'SO 1-1,3 ij 


(27) 


v. . 

i,j- 


h 




V 


ij 


) 


To achieve upstream differencing for the convection terms, the value of 
a convected property is defined as the value existing in the upstream 
cell as determined by the sign of the convecting velocity at the cell 
interface. Convected properties evaluated in this manner are p,pu,pv, 
and (pe + p) . The result is that the upstream property is convected 
across the interface at the average velocity. A typical flux at the 
left face of cell (i,j) is illustrated as 


P i-% ,j U i-*5»j 


(28) 


The difference in mass flux across cell (i,j) in the x-direction 
becomes 

6 (pu) . . = p. a. ,u. u . — P . l . , u- u i (29) 

x ij l-HSjj i+'SjJ i-1ij i- *5, J 


Since the quantities at cell walls are invarient during calculations 
at the upstream and downstream cells, local conservation principles are 
observed identically. Whatever quantity leaves one cell across a cell 
wall must enter the adjacent cell. 


An attempt was made to describe the transition of fluid properties 
from one cell to the next as a parabola which was biased in the upstream 

direction. While this scheme was numerically stable, it did not 

eliminate temperature decreases in cells surrounding the high temperature 

heater cell. 
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The pressure gradient across the node is represented by a centered 
space difference. To accomplish this in the present formulation, the 
cell wall pressure is defined as an average: 


^ (P i-l,j + P i,j } 


The resulting pressure difference at cell (i,j). 


(30) 


6 

x 



= P. . u 

1 + 



(31) 


is identical to the centered pressure difference taken across two cell 
intervals. 


2 

3 T 


3 2 u 


are not influenced by 


The diffusion terms such as - 7 — 7 - and 

3x" 

fluid convection and usual centered differences are used: 


6 2 T . 
x ij 
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i-l,j 



+ T 


i+l,j 


(32) 


While not employed here due to time limitation, it would be convenient 
to express the first temperature differences at the cell walls and 
obtain the second differences as 


6 

x 


2 



6 T j . - 6 T. , , 
x i+h,3 x 


(33) 


A temperature-dependent conductivity could be added by a simple 
extension as: 


V k T) u ’ (k 


T) i +W - (k 


6 T) . , . 

x ' i-~h > j 


(34) 


The difference forms of the governing equations employed are shown be- 
low along with the necessary supporting equations and definitions. 

Unless otherwise indicated, data is taken at the n^ time plane. 
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ij Ax x ij y ij 

pe? + ~ [-6 phu.. -6 phv +-^— ( 6 2 T + 6 2 T)3 

ij Ax x F ij y K ij Ax x y 
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Typical quantities used in the above equations are illustrated 
below in the x-direction. Quantities relating to the y-direction are 
analogous . 


ph u * + p ij 

r n+1 n+1 

x PU ij P i+ J 2,j U i+ J 2,j P i- 1 5,j U i- J 5,j 


(36) 


6 puu.. = pu, . ,u.. . - pu. , .u. . 

x ij i +h,3 1-^,3 i ~h,J 


(37) 


6 pvu.. = pv,,. .u.,, . - pv. , .u. , . 

x ij i+^.j i+h,j i~h,j 

6 phu. . = (pe , ,+P , .) u n+ } . - (pe. , .+P. , .) u n+ ,^ 

x ij i+^.J i+^.J i +h,3 i ~h,3 i~h,3 


(38) 


6 2 pu.. = 

X 1J 

- 2pu i 3 

+ pu 

6 2 pv = 

x ij 

ov t-i,j - 2pv u 

+ pv 


/ n+1 

n+1 

n+1 

_ i / PU i+l,.1 

pU i1 

u ±+h,i 


n 


(39) 


i+l,j W ij 


P i^,j = *■*« + p i+ 1#J > 


(40) 


If u 


i+ %,l 


> 0 


P i +J 5»j P ij 


If ^ ; 0 


P ±+k,3 P i+l,j 


PU i +h,3 

= PU ij 

PU. .1 . 

= PU i+l,j 

PV i-^,j 

= PV i3 

PV i+ i 5,j 

= PV i+l,J 

Pe i^,j 

= Pe i3 

Pe i^ 2 ,3 

= Pe i+l,j 


1 P 

ij (work terms) 

P i-^,j 

= P l+l,j 
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Stability Conditions 


Stability of finite-difference equations is discussed in detail 

[4] 

in Richtmyer and Morton. The basic requirements for a stable 

differencing scheme is that the difference equations converge to the 
differential equations in the limit as Ax -*■ 0, At -> 0 which implies 
that disturbances in the solution decay with time. Stability con- 
straints limit the maximum time increment (At) permissible. 

The difference analogs of the Navier-Stokes equations are too 
complex to be fully analyzed with current stability analysis methods. 
The standard procedure is to evaluate the stability of the hyperbolic 
and the parabolic terms separately and to use the more restrictive of 
the two stability constraints. 

In the explicit formulation of the Navier-Stokes equations, the 
hyperbolic limit is usually the most restrictive and results in the 
stability condition 

( | u | + c) ~ £ 1 (42) 

where u' is the fluid velocity and c is the local adiabatic speed of 
sound. If u << c, the condition may be interpreted as limiting the 
propagation of a pressure wave to the distance of one space increment 
during a time step At. 

For a two-dimensional problem, the stability condition becomes 

d u l+ c) -Jf (43) 

In the present problem, under the conditions being considered, 
the speed of sound is about 2500 ft/sec. With a grid spacing of 

-4 

Ax-.l ft. the theoretical stability limit gives At=.283 x 10 sec. 
Modification For Tank Stretch 


Tank stretch affects the fluid state primarily by changing the 
density of the fluid properties. The work done on the boundaries 
changes the internal energy and also must be included. 

During a tank expansion, the mass residing in an arbitrary volume 
becomes distributed throughout a larger volume. In the present Eulerian 
formulation, it is very difficult to adjust the x-y fluid boundaries to 
achieve this volume increase. The change was taken up in the z direction 
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by permitting the unit depth to increase. Since no fluid crosses a node 
boundary during such an expansion, the mass (m) in the volume is un- 
changed. 

2 

Multiplying through by the node volume, V=L L^, the continuity 
equation on a mass basis becomes: 


3m - 3pV _ . apu ■ 3pv 

at ~ at " ax a y ; 


( 44 ) 


Upon discretizing the time variable and introducing the forward time 
difference, the equation becomes: 


n+1 n _ n+l„n+l n„n 
m -m p V -p V 


At 


At 


-v n 

dX dy 


v n +i 

Finally, dividing by — — — and rearranging, 


n+1 V n r n . . ,3pu . 9pv N , 

° m ~j*r [p +4t ( 8r + w )] 


The above volume ratio is represented as 


,n 


,n 


str 


.n+1 


V n + AV 


( 45 ) 


( 46 ) 


( 47 ) 


where 


AV 




AP = change in average tank pressure during At. 


( 48 ) 


The same argument can be used to develop the expressions for the 
momentum and energy equations. However, the rate of work done on the 
boundary during an expansion must be included in the energy equation. 


dW 2^£ dV 
dt dt V dt 


( 49 ) 
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Upon discretization 


AW AV 

At At 


(50) 


Next, dividing by V n+ ^ as before, the work done per unit volume during 


At becomes: 

At 




AW 

PAV V n+ 1 -V n 

(51) 


V n+1 

v n+1 v n+1 


or 

AW 

v n+i 

v n 

■ p (1 - ^ > 

(52) 

and finally. 

AW 

v n+1 

- P < 1 - C str ) 

(53) 

Since this is 

the work 

done on the boundary by the fluid, 

it must be 


subtracted from the available internal energy in the cell. 


The conservation equations modified for tank stretch are collected 
together below in a form showing just the time differences: 


n+1 




pu 


n+1 


pv 


n+1 


C [pu 11 + At (. . . ) ] 

str 


C . [pv 11 + At (. . . )] 

str 


(54) 


pe n+1 = C gtr [pe n + At (...)]- P (1-C str ) 


where the average tank pressure P is used for the work term. 


For computational reasons, C gtr lags by one time step the above 
calculations so that in reality it is defined as: 

n-i 

C = - 

str „n 


(55) 
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NUMERICAL BOUNDARY CONDITIONS 


Boundary Location 

Consistent with the control- volume concept, solid boundaries were 
located at the outward-looking faces of the boundary cells. This 
location has distinct advantages in simplifying the definition and 
imposition of boundary conditions. 


The alternative to this formulation (i.e. placing boundary mesh 
points at the wall surface) creates problems in evaluating the fluid 
properties at the wall, and global conservation may not be observed 

computationally. In thio context, Roache and Mueller^ 16 -* cited the 
difficulty of imposing adiabatic (and therefore any chosen heat flux) 
conditions at the wall. They also observed that non-conservative 
methods used to define the wall density led to slow but continuous mass 
loss in a flow problem over a backstep. 

The procedure of defining flux quantities and convection velocities 
at cell faces further simplifies the imposition of boundary conditions, 
and permits a unified application of the conservation equations to both 
interior and boundary cells. This procedure also avoids the use of 
imaginary grid points located outside the fluid boundary defined by 
reflection principles. The present method is considerably simpler to 
apply when flow obstructions or wall irregularities are considered. 

Boundary Values 


Convection velocities are set to zero at boundary surfaces thus 
assuring strict global conservation for convective terms. Pressure at 
the wall, which is necessary only for the momentum equations, is 
obtained by linear extrapolation from the first two normal interior 
grid points. Thermal boundary conditions take advantage of the fact 
that the heat leak is more accurately known than the temperature of 
the inner tank wall. When calculated at a boundary cell, temperature 
difference equations are employed which make use of the relations 


3T 


q x leak 


3x 


wall 


k 


( 56 ) 


^T _ Qy le ak 

» y wall k 


where q x leak and q y leak are the boundary heat fluxes obtained by 
distributing the known heat leak rate Q leak uniformly over the exposed 
surface area. This formulation identically satisfied the global energy 



balance and avoids the difficulty in defining the thermal boundary condi- 

ri6i 

tion discussed by Roach and Mueller. 


Separate equations are used at boundary cells to evaluate the 
viscous dissipation terms in the momentum equations. For example, 
instead of the centered second difference terms used at interior cells 
which have the form: 


(S X pu) ij ■ <pu i-!,j ■ 2pu i: +pu i+i,j ) 


(57) 


forward and backward second differences are employed which make use 
of the zero velocity conditions at the wall: 

2 

(<$ x pu) ± . = (-3 pu 1> + Pu i+1 j) at left wall 

' J ‘ J lr ’ J (58) 

2 

(6 pu) = (pu - 3pu ) at right wall 


Boundary conditions are imposed at the exterior face of cell 
(20, 10) such that the prescribed mass withdrawal rate occurs. Since 
the sign of the convection velocity must be positive, the 

necessary quantities at the cell wall (P, P, pu, pv, and pe) are obtain- 
ed by linear extrapolation of the form 


p i+ ‘S.J ■>7 < - p l-l,j +3p i j ) (59) 

The necessary convection velocity is obtained from the mass flow 
relation lh=pAu as 


m 

U 204,10 “ 2 (60) 

P 204, 10 

With the cell wall quantities thus defined, the difference equations 
are applied at the outlet port cell as at any other cell. The convec- 
tion of mass»x- and y- momentum, and enthalpy from the system takes 
place automatically. 
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SCALING PRINCIPLES 


Computational stability considerations place stringent require- 
ments on the permissible size of the time step (At) required to inte- 
grate the governing equations in the explicit formulation. The time 
step is determined by the time interval required for a pressure wave 
to propagate the distance of one cell. Time steps of this order are 
required for problems in high-speed flow or for problems in which pres- 
sure wave propagation contributes significantly to the solution. In the 
present problem, wave propagation may be of importance if fluid oscil- 
lations following a G-spike are of sufficient magnitude to destroy 
thermal gradients generated by heater operation under a low-g environ- 
ment. 


However, fluid heating and stratification occur on a time scale 
very much larger than that of wave motion, and the resulting computer 
time required to simulate an appropriate duration of flight is exces- 
sive. For this reason it was necessary to introduce scaling principles 
to transform the original problem into an equivalent scaled problem in 
which certain mechanisms of importance occur in the scaled time very 
much faster than in real time. Since non-linear real gas properties 
are used to describe the fluid temperature and pressure, the thermo- 
dynamic state cannot be scaled easily. Therefore, the conditions 
governing thermal stratification and fluid mixing are adjusted so that 
these mechanisms operate in scaled time. The desired rates of heating, 
fluid withdrawal, and fluid motion are increased such that the result- 
ing pressures and temperatures remain unaltered at corresponding times 
in the scaled and unsealed systems. 

The specific heat (C ) , pressure (P) , temperature (T) , and density 
(p) are not altered. However, the transport constants, thermal con- 
ductivity (k) and absolute viscosity (y),are adjusted as required. 

Since the values of y and k for supercritical oxygen are small and do 
not dominate the stratification and mixing mechanisms, these parameters 
are considered constant. Therefore, the following principles have been 
applied. The scale factor is designated by s, and the subscript s 
represents a value appearing in the scaled system such that the real 
time t is given by 


t = st g ( 61 ) 

The following constraints have been placed upon the scaled system: 

p = P (density) 

S ' ( 62 ) 
T g = T (temperature) 
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(specific heat) 


(C) 

P s 


L 

s 



(characteristic length) 


Since it is desired that the flow processes occur more rapidly, 


u 

s 

= su 

(x velocity) 


V 

= sv 

(y velocity) 

(63) 

s 


q s 

= sq 

(heat flux) 



m = sm mass flow rate (mass flow rate) 


s 


The similarity parameters which apply are: 


Reynolds Number 


Prandtl Number 


pVL _ inertial forces 
e - y - viscous forces 

C u 

_ heat generation 
r - k - heat conduction 


Nusselt Number 


Grashof Number 


qL _ total heat transfer 
u - kAT ” conductive heat transfer 

p^ggAT _ bouyant forces 
G r " “ viscous forces 


To maintain similarity of R g ,y is scaled as: 
y g = su. 


(64) 


Similarity of both and N^ require that: 

k = sk. (65) 

s 

Recognizing that the heat diffusion rate must be scaled, the require- 
ment placed upon k can be obtained directly from the definition of the 
thermal diffusivity. 

a = sa = s ( ~ — ) . (66) 

s PC/ 
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) 


Since the compressibility (3) remains unchanged by the initial con- 
straints, to maintain Grashof number similarity, the acceleration 
must be scaled as: 



It should be emphasized that this scaling procedure speeds up the 
physical processes occurring in the fluid without disturbing the 
thermodynamic properties of the fluid itself. To accomplish this, it 
is necessary to adjust only the transport properties M and k. In par- 
ticular, the properties which determine the sonic velocity in the fluid 
have been preserved. The basic purpose for scaling is to increase the 
ratio of the fluid velocity to the sonic velocity (Mach number) 
because the sonic velocity of the explicit numerical solution limits the 
size of the permissible computation time steps. 

Verification of this scaling procedure was performed by comparing 
temperature and velocity profiles taken at corresponding times from com- 
puter runs in which differing scale factors are used. In addition, good 
results were obtained for the Apollo 12 pressure collapse and heater 
cycle simulations in which scale factors of 2400 and 6000 were employed. 

The increase in acceleration required to maintain the Grashof 
number in the scaled solution affects the required hydrostatic pres- 
sure distribution and, therefore, the distribution of the fluid mass. 
However, for problems involving a low-G environment, a scale factor 
which permits a reasonable simulation time gives rise to a negligibly 
small fictitious pressure gradient. 

It appears that the maximum valid scale factor is restricted 
by the magnitude of the acceleration. Figure 3 shows the relative 
increase in hydrostatic pressure from the center of the tank to the 
tank wall necessary to counteract the scaled gravitational body forces. 
At a given G- level, experimentation indicates that a scale factor that 
produces a 1% pressure increase does not alter the convection or thermo- 
dynamic behavior observed in the Apollo 12 simulation. It should be 
cautioned that scaling amplifies fluid-dynamic start-up transients 
resulting from imprecisely known initial conditions (eg. the velocity 
profile around the outlet port) , so that additional scaling restrictions 
must be considered. 

In addition, it should be noted that stability conditions relate 
permissible time steps to thermal conductivity and to fluid viscosity. 
However, these constraints are less restrictive of the scale factor 
than the G-constraint for the present problem. 
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NUMERICAL RESULTS 


Natural Convection in a Rectangular Enclosure 

A test case in natural convection was run and the results are 
compared with an incompressible solution of a similar problem obtained 

by Wilkes and Churchill. 

A two-dimensional rectangular enclosure containing an incompressi- 
ble fluid is oriented in a vertical plane with respect to the gravity 
vector. The left wall is held at a constant temperature T , and the 
right wall is held at a higher temperature also assumed constant. 

The other two walls are insulated. Initially, the fluid is at rest and 
at a temperature equal to the average of the boundary wall temperatures. 
The problem is to find the fluid velocity components and temperature 
at points throughout the fluid as steady-state conditions are approached. 

As heat conduction takes place, a negative horizontal density 
gradient develops. The resulting unbalanced vertical buoyant forces 
cause two vortices to form with a net counter-clockwise rotation of the 
fluid. As the flow continues, the density gradient deforms and the two 
vortices eventually merge into one due to the viscous dissipation of 
momentum. As steady-state conditions are approached, the moment pro- 
duced by the viscous forces acting at the walls balances the net buoy- 
ant moment. 

Wilkes and Churchill introduced the vorticity and stream functions 
into the non-dimensionalized equations of motion. A linear density 
dependence upon temperature was used to produce the essential density 
gradients for natural convection. An implicit alternating direction 
technique was used to advance the time-dependent solution toward steady 

state. Figures presented^^ show the stream function and isotherms 

at several times including steady state. The dimensionless conditions 

associated with these figures are: P = .733, G = 20,000, N = 2.874. 

& r r n 

The formulation used in the present solution necessitated two 
modifications to the problem. First, fluid compressibility was intro- 
duced in order to compute pressure explicity in terms of temperature 
and density. The ideal gas relation was assumed for this purpose. 
Second, the constant- temperature boundary conditions were changed to 
constant-heat-flux boundary conditions to be compatible with the model 
capability developed for simulation of heat leak into the Apollo oxygen 
tank. 
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The compressibility modification had a negligible effect on the 
solution since the vertical pressure differential was kept small with 
respect to the bulk fluid pressure. The second modification resulted 
in observable differences in the temperature profiles near the vertical 
walls. However, the basic character of the temperature profile away 
from these walls is preserved, and the resulting velocity profiles 
are similar to those obtained by Wilkes and Churchill. 

The fluid properties, heating rate, and acceleration were adjusted 
so that the proper boundary temperatures were approached at steady- 
state. The final dimensionless quantities achieved were: Pr = .611, 

= 21,000, and Nu = 3.88 which are in reasonable agreement with the 
Wilkes and Churchill values and yielded similar results. These values 
were computed using the average extrapolated wall- temperatures. 

A 10 x 10 cell grid was used to describe the fluid volume. 
Dimensions of the square enclosure, the bulk density, and the compressi- 
bility of the fluid were selected to yield a problem in which heating 
and fluid motion took place rapidly in real time in order to minimize 
computer time. Problem conditions are shown in Figure 4. 

Figure 5a is a computer-generated velocity vector plot which shows 
the two initial counter-clockwise rotating vortices after 1/2 second 
of flow development. Figures 5b and 5c show the two vortices merging 
into one and the development of circularized flow. Figure 5d shows 
the essentially steady velocities at 4.95 seconds. Figure 6 shows the 
same vector plot at 4.95 seconds on which the dimensionless stream 
function obtained from the Wilkes and Churchill solution has been super- 
imposed. Steady-state temperature profiles at ten horizontal cross- 
sections and fluid isotherms are shown in Figures 7 and 8. The general 
shape of the isotherms is in good agreement with the dimensionless 
isotherms of Wilkes and Churchill. 


It is important to note the difference in the thermal boundary 
conditions. Isotherms cannot intersect the constant temperature walls 
used by Wilkes and Churchill. However, to maintain a uniform heat flux 
at the wall as assumed in the present solution, the wall temperature 
must vary in the vertical direction and isotherms intersect these walls 
as shown. Although the present results cannot be compared directly with 
those of Wilkes and Churchill, the general character of the solutions 
appears to be in reasonable agreement. 
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It is interesting to note that starting from zero, the fluid 
rotation rate passes through a maximum before approaching the steady- 
state rate. This phenomenon is the result of the inertial lag in 
responding to a change in heat flux at the walls (as occurs at t = 0) . 
During this lag, high temperature and low temperature fluid masses 
accumulate near the respective walls, resulting in an over-acceleration 
of the fluid. The bulk rotation rate then builds up to a maximum and 
carries the heated and cooled fluid elements across the vertical center- 
line. At this point, a net counter-acting moment develops which re- 
tards the bulk rotation, and the steady-state rate is approached 
asymptotically. 


Apollo 12 Pressure Collapse Simulation 


The program was operated for flight conditions to demonstrate 
capability to simulate the stratification and mixing of supercritical 
oxygen which takes place under a flight- type acceleration environment. 
These conditions are shown in Figure 9 and are approximately those of 
the Apollo 12 mission at 7:30 GET. The density corresponds to the 95% 
tank quantity. 


The heat added to the fluid cross section was input at cell 
(12,10) as shown in Figure 9. This single heater cell represents 
a heater boundary layer volume of 1.73 cubic inches. The rate of oxygen 
withdrawal was 1.4 lbm/hr. The outlet port was the exterior face of 
cell (20,10). The pressure limits for the heater switch were set at 
860 and 900 psi so that the heater cycled automatically keeping the bulk 
pressure within a 40 psi dead band. A constant acceleration of 


”8 

2x10 G’s was applied in the -Y directio^. At a simulation time of 
70.5 minutes, an acceleration step to 10 G's was applied in the -X 
direction and the heater was turned off. The scale factor of 2400 
was used in this simulation to achieve a computer time to simulation 
time ratio of 5/1 using a program time step of .5 x 


10 ^ seconds. 


Figure 10 shows the stratified tank pressure (upper curve) and 
the equilibrium pressure (lower curve) as functions of time. The simu- 
lation started from equilibrium conditions. The equilibrium pressure 
rise rate was 3.7 psi/minute and the decay rate was -3.8 psi/minute. 
After 70 minutes of stratification, the minimum potential pressure 
collapse developed was 80 psi. The general divergence of the two curves 
shows that a quasi-steady state pressure collapse potential had not been 
reached after 70 minutes of stratification. The time required for the 
first complete heater cycle is about 12 minutes. The cycle time for 
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succeeding cycles decreased to about 3 minutes per cycle which is about 
one-third the time required for a flight heater cycle under similar condi- 
tions. The inclusion of tank stretch would about double this cycle 
time. The maximum temperature reached in this simulation was 287 °R at 
the heater node. Negligible fluid convection occurred in the vertical 
direction. 

Figure 11 shows the tank pressure decaying toward the equilibrium 
pressure in response to the G-step to 10~4 G’s at 70.5 minutes simula- 
tion time. Fluid oscillations in the tank caused by the G-step and 
accentuated by scaling have been smoothed. The pressure decay observed 
in the Apollo 12 data at 8:36 GET also is shown for comparison. 


Apollo 12 Heater Cycle Simulation 

Following the relatively successful initial simulation which 
demonstrated thermal stratification with ensuing pressure collapse, it 
was decided to incorporate refinements to the idealized model in an 
effort to more accurately predict the heater cycle times. These refine- 
ments included the effects of tank stretch, line compression, heater 
thermal mass, and heater radiation to the tank wall. 

Under similar environment and initial conditions but with somewhat 
more rigorously defined heating and mass flowrates, a second Apollo 12 
simulation was undertaken which incorporated the above refinements. 

The time step of .25 x 10“^ seconds as indicated by stability conditions 
was observed even though a time step twice as large did not appear to 
alter the numerical stability. The problem was scaled by a factor of 
6000 so that 60 minutes of simulation time was covered by 0.6 seconds 
of solution time. The associated computer time was three hours result- 
ing a computer time/real time ratio of 3/1 on the 1108 system at the 
NASA-MSC computing center. Considering the change in scale factor and 
the smaller time steps taken, revisions made to the program allow it to 
run two times faster than it did for the first Apollo 12 simulation. 

Figure 12 shows the resulting computed tank pressure historv. 
Pressure switch limits were set at 860 and 900 psi. The run was per- 
formed in a number of segments using the program re-start capability. 

At the beginning of the second and third segments at 13.5 and 23 minutes 
the heater switch was inadvertently reset to the "on" position which 
accounts for the fact that the pressure does not decay to the lower 
pressure limit at the end of the first two cycles. 
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The fact that the tank stretch and other refinements do indeed 
lengthen the heater cycle times may be seen by comparing Figures 10 
and 12. The heater cycle time appears to be around 12 minutes which 
agrees well with the 13 minute cycle time predicted from Figure 3.4.10 
of the Apollo Handbook. [1] 

The effect of heater thermal mass is observed to cause pressure 
excusions beyond the pressure limits and a rounding of the pressure 
peaks. The developed shape of the pressure decay from the peak is due 
to the decay of the local high temperature at the heater node followed 
by the more gradual decay resulting from the general fluid expansion 
caused by mass withdrawal. The highest temperature reached by the 
fluid in the heater node was 293°R. 

In this problem, fluid convection was virtually non-existent. After 
one hour of simulation, the only apparent migration of the fluid was 
toward the outlet port which is located perpendicular to the acceleration 
vector. Therefore, under 2x10“® G's these results indicate that conduction 
is the primary heat transfer mechanism. 


Acceleration Effects 


The above Apollo 12 simulation was repeated but under an acceleration 
of -2x10”® g's - two orders of magnitude higher - to illustrate the 
effects of acceleration on convection velocities. Employing a scale 
factor of 2400, the simulation was run to 11.5 minutes by which time a 
maximum convection velocity of .35x10”® ft/sec had developed at the heater 
node. Figure 13 shows a velocity vector plot of the convection pattern 
developed by time. Density and temperature data indicate that the 
heated fluid migrated a distance of about 1-1/2 cells or 1.8 inches 
which corresponds to an average velocity of about .2x10”® ft/sec. 


Scaling Verification 

It was tacitly assumed in the discussion of Scaling procedures that 
all the transport process would in fact take place in scaled time. Also 
assumed was that the pressure and temperature distributions would not be 
altered by the scaling procedure. 

The Apollo 12 convection problem just described was repeated using a 
scale factor of 4800 instead of 2400. Comparisons of the data from the 
two cases indicate that the temperature profiles and pressure rise 
rates are preserved. At a scaled time of 5 minutes for both cases the 
tank pressures agreed within 1 psi and the heater cell temperatures agree 
within 1°R after rising 53°R during heater operation. Scaled fluid 
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velocities in the second case (scale factor = 4800) are just double the 
scale velocities of the first case (scale factor = 2400) as required. 
Upon descaling, the actual fluid velocities obtained from the two scaled 
systems are the same, and the scaling procedure is verified. 

Evidence supporting the scaling procedure also can be inferred 
from the quality of the Apollo 12 simulations which were scaled between 
3 and 4 orders of magnitude. 


High Heat Leak and Boundary Roughness 

The circular cross sectional geometry of the Apollo oxygen tank 
was approximated in a step-wise fashion. A test case was performed to 
investigate the affect of such a boundary on the flow pattern. In this 
case, high heat leak of 2000 BTU/hr was imposed at the fluid boundary. 
This rate corresponds to the limiting heat leak that would occur if the 
annulus vacuum were lost. [3] An acceleration level of 10 -5 G's was 
imposed in the -y direction, and the scale factor of 4800 was used. 

The velocity vector plots in Figure 14 show the natural convection 
^fter 10, 20, and 32.5 minutes of simulation. The maximum velocity 
indicated is 10“ 5 ft/sec. Figure 14c shows that local flow distortions 
are introduced at the protruding corners. This surface roughness 
probably makes this approach unsuitable for investigations involving 
tank rotation. 


SUMMARY 


The conservation equations governing the motion of a compressible 
viscous fluid were solved in two dimensions using an explicit finite- 
difference technique. The difference equations were formulated in terms 
of control- volume grid cells which simplified the imposition of heat 
flux boundary conditions and assured computational observance of local 
and global conservation principles. This system also permitted the 
unified application of the difference equations to both interior and 
boundary cells without resorting to exterior cells and reflection 
principles. Real-fluid properties describing the thermodynamic behavior 
of supercritical oxygen were used so that the pressure collapse phenomenon 
could be observed in the Apollo oxygen Cryogenic Storage System operating 
under low - G conditions. 

The numerical procedure was applied to the simulation of thermal 
stratification and fluid mixing in the Apollo oxygen storage tank. 
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The convex tank geometry was approximated by removing the corner cells 
of a rectangular cell grid* Wall boundary conditions were adjusted 
at one cell face to emit a prescribed mass flowrate. Electrical heater 
input was treated as local internal heat generation. The effects of 
heater thermal mass and radiation to the tank were included. The effect 
of tank stretch and line compression on dP/dt was modeled as an out-of- 
plane fluid expansion. Scaling principles were invoked to achieve 
acceptable computer execution times for reasonable flight durations. 

A verification test case was performed involving heat transfer and 
natural convection in a vertically-oriented rectangular volume. The 
convection pattern and isotherms are in good agreement with another 
numerical solution. 

A simulation of stratification and mixing occurring around 7:30 GET 
in the Apollo 12 mission was presented. Natural convection under 2 x 10 
G's acceleration was shown to be negligible. The pressure collapse of 
about 75 psi following a simulated vehicle maneuver is compared with 
flight data with good results. 

Modifications in the heater characterization along with tank stretch 
effects were shown to significantly improve the simulation of heater 
cycle operation. A number of additional cases were presented to show 
the effects of higher acceleration on convection velocities to verify 
the scaling techniques employed and to evaluate the effects of boundary 
roughness on convection patterns. 

The results obtained for these initial test cases indicate the 
general capability of this analysis. Unfortunately, time limitations 
prevented the refinement of certain empirical considerations which would 
further improve the accuracy of the Apollo simulations. 

Future Developments 


It has been observed that the finite-difference solution of partial 
differential equations is limited more by currently available theoreti- 
cal understanding than by computer capability . ^ J However present 
computers, which perform all operations serially, are not particularly 
well adapted to computing finite-difference solutions which proceed in 
a series— parallel fashion. It seems reasonable that the multi— process- 
ing capability of present machines will be extended to parallel process- 
ing within the same computer program. Such is the thrust of the new 
multi-computing system at the University of Illinois. ^ This system, 

called ILL I AC IV, consists of 64 independent processors which operate in 
unison. By using parallel processing, computer run time can be reduced 
by a factor of 64. Thus it appears certain that next generation com- 
puters of this type will significantly increase the capability of finite 
difference solutions in multiple space dimensions. 
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NOMENCLATURE 


i 


Symbol 

C. 

lag 

C rad 

C str 

C . 
vol 

e 

E 

g 

k 

P 

q 

Q 

^leak 


v 

t, At 


u 

v 

V 

P 

P 

pe 

pu 

Pv 

Ax, Ay, L 

T 

Subscript 

i 

j 

x 


y 

Superscript 

n 


Description 

heater on/off ramps (0 _< £ 1) 

heater radiation factor (0 ^ ^ ra( j — 
tank stretch factor ( = V n ^/V n ) 
model/actual volume ratio (Apollo 12 = .075) 
internal energy per unit mass 

2 2 

total energy per unit mass ( = e + % (u + v )) 

gravitational acceleration 

thermal conductivity 

pressure 

heat flux 

internal heat generation rate 
boundary heat leak rate 

internal heat generation rate per unit volume 

time, time increment 

temperature 

velocity, x-component 

velocity, y-component 

tank volume 

absolute viscosity 

density 

internal energy per unit volume 
momentum, x-component 
momentum, y-component 
node dimensions 
shear stress 

x-direction index 
y-direction index 
x-component 
y-component 

time step index 
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Figure 1 Oxygen Temperature - Energy Relation 
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Direction 


Figure 2 Grid System and Indexing Procedure 
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i-l» 1+1 

• 

i, j+1 

• 

i+1, j+1 

• 
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• 
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1 > 1 

- • i 

i-% ,1 1+ 

1,1-H 

i+1, j 
t, • 

i-1. .1-1 

• 

i> 1-1 

• 

i+1, 1-1 

• 


x - Direction 


• - location of cell properties 
X - location of fluxes 
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SCALE FACTOR 



Figure 3 Scale Factor vs Acceleration 






Figure 4 Natural Convection Conditions 



P q =49.71 lbm/ft^ 
z = .3 

MW = 32. lbm/lbm-mole 

C = .2 B/lbm-°R 
v 

y = .45 Poise 
k = 39 B/f t-hr-°R 
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Figure 7 Steady-State Temperature Profiles 



Figure 8 Steady-State Isotherms 


ISOTHERM PLOT AT T = 7.5 SECONDS 
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Figure 9 Apollo 12 Simulation Conditions 
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ACCELERATION, a 

Conditions 

P Q = 867.7 psi Heater = 1 node 

T q = 200 °R q = 528.6x(l/17) B/hr 

p Q = 65.45 lbm/ft 3 w = 1.4x(3/40) lbm/hr 

= -2. x 10" 8 G 's Scale = 2400 

17U 


a 


PRESSURE (PSI) 



Figure 10 Apollo 12 Stratified and Equilibrium Pressures 










PRESSURE (PSI) 



Figure 12 Apollo 12 Stratified Pressure. - Including Refinements 
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Figure 14 Surface Effects Under High Boundary Heat 
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